Chloe Bracis
What is the animal doing? Tools for exploring behavioural structure in animal movements
Eliezer Gurarie, Chloe Bracis, Maria Delgado, Trevor D. Meckley, Ilpo Kojola and C. Michael Wagner
Journal of Animal Ecology, 2015
Questions
Tools
Application
Explotatory
Explanatory
Predictive
Supervised
Unsupervised
Make use of information you have about the system!
require(waddle)
data(Lamprey)
# smooth track for GPS error when stationary
Lamprey.smooth = SmoothTrack(Lamprey, 3)
FPT - time it takes for an individual to enter and leave a circle of fixed radius r drawn around each location
Parameter to specify: radius
require(adehabitatLT)
Lamprey.traj = as.ltraj(data.frame(Lamprey$X, Lamprey$Y), Lamprey$Time,
id = "Lamprey")
Lamprey.smooth.traj = as.ltraj(data.frame(Lamprey.smooth$X, Lamprey.smooth$Y), Lamprey.smooth$Time,
id = "Lamprey")
radii = c(5, 10, 20)
Lamprey.fpt = fpt(Lamprey.traj, radii)
Lamprey.smooth.fpt = fpt(Lamprey.smooth.traj, radii)
Other implementations:
modpartltraj in adehabitatLT and (Guéguen 2001, 2009)Assumptions:
First, regularize data:
Lamprey.smooth.reg = InterpolatePoints(Lamprey.smooth, 2, "min")$Data
Lamprey.smooth.traj = as.ltraj(data.frame(Lamprey.smooth.reg$X, Lamprey.smooth.reg$Y),
Lamprey.smooth.reg$Time, id = "Lamprey")
See modpartltraj in adehabitatLT
Lamprey.smooth.segments = Prep.segments(Lamprey.smooth.traj, units = "min",
sd = 5, nmodels = 20)
Maximum likelihood for K = 4
See partmod.ltraj in adehabitatLT
Lamprey.smooth.partition = Partition.segments(Lamprey.smooth.segments)
plot.segments(Lamprey.smooth.partition,
xlab="", ylab="Step length (m)")
DiagPlot.segments(Lamprey.smooth.partition)
Smooth
Flat
Format data and set tuning parameters
Lamprey.VT = GetVT(Lamprey.smooth, units = "min")
windowsize = 30
windowstep = 1
K = 0.5
Choose function of data columns for analysis
Lamprey.ws = WindowSweep(Lamprey.VT, "V*cos(Theta)",
windowsize, windowstep,
plotme = FALSE, K = K,
tau = TRUE, progress = FALSE)
plot(Lamprey.ws, type="smooth", threshold = 2,
legend = FALSE)
plot(Lamprey.ws, type="flat", clusterwidth = 3,
legend = FALSE)
DiagPlot(Lamprey.ws, "smooth")
DiagPlot(Lamprey.ws, "flat")
Movement model: correlated random walk (CRW)
Assumptions:
Model fitting:
Movement:
State switching:
require(moveHMM)
Lamprey.hmm = prepData(Lamprey.smooth, type = "UTM", coordNames = c("X", "Y"))
plot(Lamprey.hmm, compact = TRUE)
summary(Lamprey.hmm)
Movement data for 1 animal:
Animal1 -- 430 observations
Covariate(s):
Time
Min. 25%
"2010-05-02 14:50:40 CEST" "2010-05-02 17:32:15 CEST"
Median Mean
"2010-05-02 20:13:20 CEST" "2010-05-02 20:22:12 CEST"
75% Max.
"2010-05-02 23:15:30 CEST" "2010-05-03 02:05:20 CEST"
Z
Min. 25% Median Mean
728016+5043025i 728631+5042149i 728632+5042149i 728671+5042191i
75% Max.
728639+5042148i 729223+5041914i
Depth
Min. 25% Median Mean 75% Max.
0.2926667 5.0550000 5.0550000 4.4891907 5.0550000 5.7146667
summary(Lamprey.hmm)
Movement data for 1 animal:
Animal1 -- 430 observations
Covariate(s):
Time
Min. 25%
"2010-05-02 14:50:40 CEST" "2010-05-02 17:32:15 CEST"
Median Mean
"2010-05-02 20:13:20 CEST" "2010-05-02 20:22:12 CEST"
75% Max.
"2010-05-02 23:15:30 CEST" "2010-05-03 02:05:20 CEST"
Z
Min. 25% Median Mean
728016+5043025i 728631+5042149i 728632+5042149i 728671+5042191i
75% Max.
728639+5042148i 729223+5041914i
Depth
Min. 25% Median Mean 75% Max.
0.2926667 5.0550000 5.0550000 4.4891907 5.0550000 5.7146667
weibullShape = c(3, 0.5)
weibullScale = c(10, 10)
wcauchyMu = c(0, pi)
wcauchyRho = c(0.7, 0.2)
doubleStateModel = fitHMM(Lamprey.hmm, nbStates = 2,
stepPar0 = c(weibullShape, weibullScale),
anglePar0 = c(wcauchyMu, wcauchyRho),
stepDist = "weibull", angleDist = "wrpcauchy")
Value of the maximum log-likelihood: -1397.219
Step length parameters:
----------------------
state 1 state 2
shape 2.5679 0.9615813
scale 38.7629 1.0270557
Turning angle parameters:
------------------------
state 1 state 2
mean -0.01419219 -0.2208763
concentration 0.85071975 0.3227927
Regression coeffs for the transition probabilities:
--------------------------------------------------
1 -> 2 2 -> 1
intercept -4.035819 -5.045129
Transition probability matrix:
-----------------------------
[,1] [,2]
[1,] 0.982635653 0.01736435
[2,] 0.006399412 0.99360059
Initial distribution:
--------------------
[1] 9.999994e-01 6.328919e-07
$stepPar
$stepPar$lower
state 1 state 2
shape 2.252089 0.8857300
scale 35.954625 0.9085821
$stepPar$upper
state 1 state 2
shape 2.927998 1.043928
scale 41.790508 1.160977
$anglePar
$anglePar$lower
state 1 state 2
mean -0.06121604 -0.4458002
concentration 0.81355297 0.2507313
$anglePar$upper
state 1 state 2
mean 0.0317344 -0.002439682
concentration 0.8805810 0.390922736
$beta
$beta$lower
1 -> 2 2 -> 1
intercept -5.433848 -6.435466
$beta$upper
1 -> 2 2 -> 1
intercept -2.637791 -3.654792
Decoding states sequence... DONE
Decoding states sequence... DONE
Computing states probabilities... DONE
weibullShape = c(3, 1, 1)
weibullScale = c(20, 1, 2)
wcauchyMu = c(0, pi, 0)
wcauchyRho = c(0.9, 0.2, 0.6)
tripleStateModel = fitHMM(Lamprey.hmm, nbStates = 3,
stepPar0 = c(weibullShape, weibullScale),
anglePar0 = c(wcauchyMu, wcauchyRho),
stepDist = "weibull", angleDist = "wrpcauchy")
Decoding states sequence... DONE
Decoding states sequence... DONE
Computing states probabilities... DONE
for (i in 1:n)
{
weibullShape = runif(nstates, 0, 50)
weibullScale = runif(nstates, 0, 50)
wcauchyMu = runif(nstates, -pi, pi)
wcauchyRho = runif(nstates, 0, 1)
startVals[i,] = c(weibullShape, weibullScale, wcauchyMu, wcauchyRho)
tryCatch(
{
model = fitHMM(Lamprey.hmm, nbStates = nstates, stepPar0 = c(weibullShape, weibullScale),
anglePar0 = c(wcauchyMu, wcauchyRho),
stepDist = "weibull", angleDist = "wrpcauchy")
angleMean[i,] = model$mle$anglePar[1,]
angleCon[i,] = model$mle$anglePar[2,]
stepShape[i,] = model$mle$stepPar[1,]
stepScale[i,] = model$mle$stepPar[2,]
nll[i] = model$mod$minimum
},
error = function(e) print(paste("i =", i, conditionMessage(e)))
)
}
AIC to choose among models
AIC(doubleStateModel, tripleStateModel)
Model AIC
1 tripleStateModel 2783.884
2 doubleStateModel 2816.437
However, AIC will frequently select more states than may be biologically relevant or interpretable.
Clustering method (default is speed and turning angle, but can specify other variables) using Gaussian mixture model (Garriga et al. 2016)
require(EMbC)
# here we use the speeds and turning angles we already calculated,
# you can also pass in a move object directly
clustering = embc(as.matrix(Lamprey.VT[, c("V", "Theta")]))
... computing starting delimiters
[1] 0 -0.0000e+00 4 428
[1] 1 -6.1054e+00 4 320
[1] 2 -5.9219e+00 4 96
[1] 3 -5.2132e+00 4 41
[1] 4 -4.3543e+00 4 39
[1] 5 -4.1245e+00 4 24
[1] 6 -4.0369e+00 4 12
[1] 7 -4.0044e+00 4 13
[1] 8 -3.9879e+00 4 5
[1] 9 -3.9819e+00 4 6
[1] 10 -3.9785e+00 4 2
[1] 11 -3.9717e+00 4 7
[1] 12 -3.9623e+00 4 6
[1] 13 -3.9574e+00 4 1
[1] 14 -3.9559e+00 4 1
[1] 15 -3.9550e+00 4 2
[1] 16 -3.9549e+00 4 4
[1] 17 -3.9544e+00 4 2
[1] 18 -3.9541e+00 4 0
[1] 19 -3.9538e+00 4 1
[1] 20 -3.9537e+00 4 0
[1] 21 -3.9536e+00 4 2
[1] 22 -3.9535e+00 4 0
[1] 23 -3.9535e+00 4 0
[1] ... Stable clustering
Nu.sim, Tau.sim, BCRW.sim in waddle, apply methodsdata(Multipaths)
# ltraj object
Nu.traj = as.ltraj(data.frame(X = Re(Nu.sim$Z),
Y = Im(Nu.sim$Z)),
Sys.time() + 1:length(Nu.sim$Z),
id = "Nu.sim")
# bcpa object
Nu.VT = GetVT(data.frame(Z = Nu.sim$Z,
Time = Sys.time() + 1:length(Nu.sim$Z)),
units = "sec")
# moveHMM object
Nu.hmm = prepData(Nu.traj[[1]], type = "UTM", coordNames = c("x", "y"))
For this section, use your own data, the Wolf data in the waddle package or another dataset
Daily positions
data(Wolf)
Wolf2 = data.frame(X = tapply(Wolf$X,
substr(Wolf$Time, 1, 10), mean),
Y = tapply(Wolf$Y,
substr(Wolf$Time, 1, 10), mean))
Wolf2$Time = as.POSIXct(row.names(Wolf2))
Wolf.traj = as.ltraj(data.frame(Wolf2$X, Wolf2$Y),
as.POSIXct(Wolf2$Time),
id = "Wolf")